arXiv:1505.06176vl [math.OC] 20 May 2015 


Numerical testing in determination of sound 
speed from a part of boundary by the 

BC-method 

M.I.Belishev* I.B.Ivanovj I.V.Kubyshkin| V.S.Semenov§ 


Abstract 

We present the results of numerical testing on determination of the 
sound speed c in the acoustic equation uu — c^Au = 0 by the boundary 
control method. The inverse data is a response operator (a hyperbolic 
Dirichlet-to-Neumann map) given on controls, which are supported on 
a part of the boundary. The speed is determined in the subdomain 
covered by acoustic rays, which are emanated from the points of this 
part orthogonally to the boundary. The determination is time-optimal: 
the longer is the observation time, the larger is the subdomain, in 
which c is recovered. The numerical results are preceded with brief 
exposition of the relevant variant of the BC-method. 
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Introduction 


1.1 About the method 

The boundary control method (BCM) is an approach to inverse problems 
based on their relations with control and system theory 13 m [n]. It is 
a rigorously justihed mathematical method of synthetic character: Rieman- 
nian geometry, asymptotic methods in PDE, functional analysis and operator 
theory are in the use. Beginning on its foundation in 1986 [1], there was a 
question whether such a purely theoretical method is available for numerical 
implementation. The first affirmative results were obtained by V.B.Filippov 
in two-dimensional problem of the density p = determination via the 
spectral inverse data [12] . Later on, an algorithm based on the spectral vari¬ 
ant of the BCM was elaborated and tested by S.A.Ivanov and V.Yu.Gotlib 
in [m |7|. 

A dynamical variant of the BCM deals with time-domain inverse data that 
is a response operator (hyperbolic Dirichlet-to-Neumann map). It provides 
time-optimal reconstruction: the longer is the observation time, the bigger 
is the subdomain, in which the parameters are recovered. It is the feature, 
which makes this variant most relevant for possible applications to acoustics 
and geophysics. The corresponding algorithm was elaborated and tested by 
V.Yu.Gotlib in [10]. It recovers the density in a near-boundary layer from 
the data given on the whole boundary. 

Time-optimal determination of density via the spectral and time-domain 
inverse data given on a part of boundary is proposed in [5]. The procedure 
uses singular harmonic functions; its spectral variant was realized numerically 
(see [5], section 7.7). In [6] and [8], its dynamical variant was modihed to 
make it more prospective for applications in geophysics, the modihcation 
being based on geometrical optics. 

In beginning of 2000’s, L.Pestov proposed a version of the BCM, which 
determines some intrinsic bilinear forms containing parameters under recon¬ 
struction via the inverse data and, then, recovers the parameters from the 
forms. This version is not time-optimal but, on expense of big enough ob¬ 
servation time, provides more stable numerical algorithms. The results of 
the collaboration, which develops this approach in the I.Kant Baltic Federal 
University (Kaliningrad, Russia), are presented in [T9|[20l|2T]. 

Recently, L.Oksanen applied the BCM for numerical reconstruction of 
the obstacle [TH] . 


2 


There also exists a time optimal and data optimal approach by V.Romanov 
[22] bnt it is not implemented and tested yet. Another (not optimal) direct 
reconstrnction methods, which are nnmerically (and experimentally) tested, 
see in 


1.2 Inverse problem 

The goal of onr work is to elaborate the BC-algorithm for time-optimal de¬ 
termination of the sonnd speed via the time-domain inverse data given at a 
part of bonndary, and test it in nnmerical experiment. 

• Let C M” be a (possibly, nnbonnded) domain with the bonndary T. 
We deal with a dynamical system 


Utt — c^Am = 0 

in 12 X (0,T) 

(1.1) 

u\t=Q = Ut\t=Q = 0 

in 12 

(1.2) 

u = f 

on F X [0, T], 

(1.3) 


where c = c{x) is a smooth enongh positive fnnction {speed of sound), f is 
a boundary control, u = u^{x,t) is a solntion {wave). With the system one 
associates a response operator 




f "^^Irxio.T] 


(1.4) 


where {■■■)u is a derivative with respect to the ontward normal z/ on T. In a 
general form, the inverse problem is to answer the qnestion: To what extent 
does the response operator determine the sonnd speed into the domain? Also, 
the determination procednres are of principal interest. 


System (1.1)-(1.3) is hyperbolic and, as snch, obeys the finiteness of the 
domains of influence (FDI). It describes the waves propagating with hnite 
speed c, and the relevant setnp of the inverse problem mnst take this property 
into acconnt. Snch a setnp is given below, after geometric preliminaries. 

• The sonnd speed indnces a travel time metric = c~‘^\dx\^ (shortly, 
c-metric) and the corresponding distance T{x,y) in fl. For a snbset A C 12, 
by _ 

12^ := {x G 12 I t{x. A ) < ^} 

we denote its c-metric neighborhood of radius f. 
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By we denote a geodesic in c-metric (ray), which is emanated from 
7 G r into in direction —u, and is of the c-length Let a C L be a part 
of the boundary. A set 

■yGcr 

is called a ray tube. On Fig la,b, the neighborhood and tube are 


contoured by the closed lines {1,2, 
(Sj is shaded). 



, 4, 5,6,1} and {5, 6, 2, 3,5} respectively 





Figure 1; Tube and domain D'^ 

If T is small enough then the ray held is regular in the tube. Let be 
the inhmum of T’s, for which such a regularity does occur. 

Convention 1. In what follows, unless otherwise specified, we assume that 
a is diffeomorphic to a disk {p G | \p\ < 1} and T <T„. Such a case is 
said to he regular. 

The part a determines the space-time domains 

:= {{x,t)\x G Oj, 0 < t{x, a) < 2T — t} and 
El := {(a;,t)|x G 0 < t < r(a:,cr)}, 

and the space-time surfaces 

0^= n (F X [0, T]} , Of n (F X [0, 2T]}. 

All of them are mapped by the projection (x, t) i—)■ a; to klf. Domain D'^ 
is shown on Fig l.b (shaded). Domain Ef. C lies under the surface 
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{{x,t) I t = r(x,cT)}, which consists of three parts countered by the closed 
lines { 6 , 7, 8 , 6 }, {5, 6 , 8 , 9, 5}, and {5, 9,10, 5}. The surfaces 0^^ and 0^ are 
countered by the lines {1, 6 , 5,4,10,13,14, 7,1} and {1, 6 ,5,4,10,11,12, 7,1} 
respectively. 

If c < c* = const holds in then for the sets 
cr^ := {7 G r I r( 7 , a) < and af := {7 G T | dist]Rn(7, a) < (1.5) 

one has <Z af, and the relations 

el = a^x[0,T] Cafx[0,T] (1.6) 


are valid. 

• Assign a control / to a class if supp/ C a x [0,2T], i.e., it acts 
from a during the time interval 0 < t < 2T. Owing to the FDI, an extension 


of system (1.1)-(1.3) of the form 


Utt — = 0 

M = 0 

u = fe:Ff 


in 

in El 


2T 


(1.7) 

(1.8) 

(1.9) 


turns out to be a well-posed problem, its solution being determined by 
the values of the speed c in the subdomain (does not depend on cIq^qt). 
The same is valid for the response operator 

( 1 - 10 ) 

associated with this problem: it is also determined by cI^t. 

By the latter, the relevant setup of the inverse problem is: for a fixed 
T > 0, given the operator Rl^ determine the speed c in flf. 

The use of the doubled time 2T is quite natural by kinematic reasons. 
The subdomain Qf is prospected with waves initiated at a. To search the 
whole Qf, the waves have to £11 it (that takes T time units) and return 
back to the boundary (for the same time T) to be detected by the external 
observer, which implements measurements at T. 

Convention 2. The operator R?J" is introduced so that, for the times 0 < 
t < T the images f are defined on the set 0^ only. For convenience of 
further formulations, we put to he extended from Qf toT x [0,T] 

by zero. 
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1.3 Results and comments 


• Let cr C r and T > 0 be given. Our a priori assumptions are that T < T^- 
(i.e., we deal with the regular case) and the sound speed upper bound c* is 
known. Under these assumptions, we propose a procedure, which recovers 
the speed c in the tube via the operator R^. Then, we demonstrate the 
results of numerical testing of the algorithm based on this procedure. 

• In fact, the procedure utilizes not the complete operator R^ but some 
information, which it determines. Namely, as will be seen, to recover c|^y, 
it suffices for the external observer to possess the following options: 

1. for any f,gG obeying the oddness condition 

/(•, T) = -/(•, 2T - f), T) = -g{-, 2T - t), 0 < f < 2T, 

one can compute the integral 

[ u-^^{7,'t)9hR)dTdt= {Rl^Jf,g)^2T, ( 1 . 11 ) 

Jax[0,2T] 

where J : —?■ R'^ is an integration: := f{-, s) ds. 

2 . for any odd / G R^'^, one can detect i-®-’ 

implement the measurements on aj (but not on the whole T!) during 
the time interval [0,T] (but not [0,2T]!) 


• In principle, the proposed procedure is identical to the versions j6| and 
[ 8 ]. Therefore, its exposition is short: we omit some proofs and derivations, 
referring the reader to the mentioned papers for detail. In the mean time, 
here we deal with more rehned (rigorously time-optimal) data that is the 
operator R"^, in contrast to [ 6 ] and [ 8 ], where the operator corresponding 
to system (1.1)-(1.3) with the hnal time t = 2T, is used as the inverse data. 


• One of the features and advantages of the BCM is that it reduces non¬ 
linear inverse problems to linear ones. In particular, the main fragment of the 
algorithm, which recovers c, is the solving a big-size linear algebraic system. 
The matrix of the system is of the form ^ ^i^h enough 

set of controls /j. As a consequence of the strong ill-posedness of the above 
stated inverse problem, this system also turns out to be ill posed but the 
linearity enables one to apply standard regularization devices. 
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2 Geometry 

2.1 c-metric 

• Let f2 C {n > 2) be a domain with the C^-smooth boundary L. A 
sound speed is a function c G C^(r2) provided c > 0. If is unbounded, we 
assume c < c* = const. 

The sound speed determines a c-metrie in hi with the length element 
= c~‘^dl‘^ and the distance 

t{x, y) := inf 

where dl is the Euclidean length element, and the inhmum is taken over 
smooth curves, which he in hi and connect x with y. In dynamics, the value 
t{x, y) is a travel time needed for a wave initiated at x to reach y. 

• Let a C L; a function 

Ta{x) := inf T{y,x ), x E Q 

yea 

is called an eikonal. Its value is the travel time from a to x. A set 

:= {x e I T„{x) < ^} (^ > 0) 

is a c-metric neighborhood of a of radius In dynamics, the waves initiated 
on a at the moment f = 0, £11 up the subdomain at t = The filled 
domains are bounded by the eikonal level sets 

r| := {x e I T„{x) = O 

(the surfaces c-equidistant to a: see the dotted line on Fig la), which play 
the role of the forward fronts of waves propagating from a into hi. 

2.2 Ray coordinates 

• Fix a point 7 G a. Let x(^,^) G be the endpoint of the c-metric geodesic 
(ray) starting from 7 orthogonally to T and parametrized by its c-length 

Also, put x(^, 0 ) = 7 . 
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For T > 0, the rays starting from a, cover a tube 

= U a;(7,0 C 

7So- (7,^)ecrx [0,T] 

In the regular case, is diffeomorphic to the set 

:= <tx | 0 ,r] 

via the map 3 (75 0 G -Bj (on Fig l.b, Sj is countered 

by {6,5,11,12,5} ). This enables one to regard a pair as the ray 

coordinates of the point G -Bj. 

• Let TT* be the Cartesian coordinate functions: 7r^{x) := x* for x = 
G The map 

(7,0^KW7,0)}r=i, (7,0 (2.1) 

realizes the passage from the ray coordinates to Cartesian ones. 

Fix a 7 G a. The equality 


c(x(7,0) = ■ 

represents c on the ray r 
in the whole tube Sj. 


E 


2 = 1 




21 2 


0 < ^ < T (2.2) 


Varying 7 , we get the sound speed representation 


2.3 Images 

• Fix a point 7 G a, choose a small e > 0, and dehne the surfaces 

^^07, 0 := {^(Y, 0 e I r( 7 ', 7 ) < £> , 0 < ^ < T. 


A function 


</(7,0 := lim 


1^^07,01 


^^0 |c^£(7,0)| ’ 


( 7,0 e so 


where |...| is a surface area in is said to be a ray spreading at the point 

a^(7,0- 



In the regular case, the coefficients 

>^(7,0 ■= := W 7 , 0 ) x ( 7 , O ]^ , 

which enter in the well-known geometrical optics relations (see, e.g., mm), 
are the smooth functions on 

• Let y he & function on ; a function y of the form 

^( 7,0 := 13(7,0 yixi7,0), ( 7,0 e Sj 


is called an image of y. For the function r(x) = 1, one has = /3. 

In terms of images, relations (2.1) and ( 2.2[ ) take the form of the repre¬ 
sentations 


(7.?)es7 


E 

2=1 


d 


^°( 7,0 


-\ 2 


(2.3) 


which will be used for determination of c in the inverse problem. 


3 Dynamics 

In section the regularity condition T < To- is cancelled, and T > 0 is 
arbitrary. However, for the sake of simplicity, we keep a to be diffeomorphic 
to a disk. All the functions, spaces, operators, etc are real. We denote 
S® := a X [0,s]. 


3.1 Spaces and operators 


Denote the dynamical system associated with problem (1.1)-(1.3) by . In 


what follows, we deal with its subsystem corresponding to controls acting 
from a. We consider it as a separate system, denote by aj, and endow with 
standard control theory attributes: spaces and operators. All of them are 
determined by c|^j,. 

• The space of boundary controls Tj := T 2 (Sj) with the inner product 


{f,9)jrT-.= f{^,t)g{^,t)drdt 
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(fir is the Euclidean surface element on the boundary) is called an outer 
space of system a^. It contains an increasing family of subspaces 

:= {/ G J-J I supp/ C a X [T - e,T]} , 0 < ^ < T 

= {0}, formed by the delayed controls acting from a. 

Here, T — ^ is the value of delay, ^ is an action time. 

• The space "Hj := c~‘^dx) with the inner product 

{y,w)^T-.= [ y{x)w{x)-^ 

JnT c^[x) 

is said to be an inner space of the system. It contains a family of subspaces 
Vi := ^y eVl\ supp?/ C , 0 < ^ < T 


{Vi := {0}), which increase as a extends and/or ^ grows. 

• In the system aj, an ‘input—)-state’ correspondence is described by a 
control operator : Vi —?• "Hj, 


W^f := u^(-,T), 


where is a solution to (1.1)-(1.3). Operator W'^ is bounded [5]. 

Since the waves governed by the equation (1.1) propagate with the hnite 
speed c, for controls acting from a one has 


suppM-^(-, .^) C , 0<.^<T. 


(3.1) 


As is easy to recognize, (3.1) is equivalent to the embedding 

W^Vi^GVi, 0<e<T. (3.2) 

• Recall that v is the outward normal to T, and the sets are dehned 


in (1.5). Denote := T x [0,T] 

An ‘input—^-output’ correspondence is realized by the response operator 
: Vi ^ L2iJ:^;dTdt), 

■= 

dehned on the set DomR^ = {/ e | /|g^x[o,r] = •^L=o = 

where is the Sobolev class and da is the boundary of cr in T). Relation 

(3.1) implies 

suppu^ C {(7,0 I 7 e cr^, 0 < ^ < T} C X [0 ,T]. 
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By the latter, for controls / G one has 


suppi?^/ (Z X [0,T] c (jf X [0,T]. (3.3) 

One more (extended) response operator —)■ dTdt) is 


f ■= 


I/|©2T , 


on 

L2 


where is a solution to extended problem (1.7)-(1.9). It is dehnec 
Domi?f = {/ e I /|g^,^[ 02 r] = /lt=o = 

is determined by the values of the sound speed c in the subdomain 
Therefore, it is reasonable to regard it as an intrinsic obiect of system aj 
(but not af!). 

Let the controls / G Domi?^ in (1.3) and / G DomiT^in (1.9)_be such 
that / = 

also coincide for the same times: 


Then, the solutions to problems (1.1)-(1.3) and (1.7)-(1.9 


y^f — yf ^ j-g^ ^ 

As a consequence, passing to the normal derivatives on T, one gets 

R^f = RYf ons^nef. 

• A connecting operator of the system is C'^ : —)■ 


(3.4) 


C‘ := {W'^yw 


The dehnition implies 

(n^(-,T),n^(-,T))^, = {W^f,W^g)yr = {C^f.g)^ , (3-5) 

i.e., C'^ connects the Hilbert metrics of the outer and inner spaces. 

A signihcant fact is that the connecting operator is determined by the 
response operator in a simple explicit way. Namely, the representation 

= 2-\S^yR^/JS^ (3.6) 

is valid, where the map S'^ : — )■ extends the controls from a x |o,r] 

to (j X [0, 2T] by oddness with respect to t = T: 


{S^f) i;t) 


f{;t), 0<t<T 

-/(•,2T-t), T<t<2T 


(3.7) 
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and J : —)■ is an integration: (Jf)(-,t) = f^f(-,s)ds (see E]-[S]). 

Note that / G BomR^ implies f G BomRf. 

As a conseqnence, we get 

(Cp.s)^ ® 'S' 2-‘ (RfjS^f,S'^g)^„ = 

0'2-i/J[Sp,S^9] (3.8) 

for arbitrary / G Domi?^ and g G 


3.2 Wave bases 


For the BCM, the fact of crucial character is that the embedding (3.2) is 
dense: the equality 


clos = Hi, 0 < e < T (3.9) 


(the closure in H^) is valid and interpreted as a local boundary controllability 
of system (!.!)-(1.3). It shows that the waves constitute rich enough sets 
in the subdomains which they £11 up. In particular, by this property, any 
square-summable function supported in can be approximated (with any 
precision) by a wave owing to proper choice of the control / acting 

from a piziEiEn]. 

• An important consequence of controllability is existence of wave bases. 

Fix a ^ G (0,T]. Let a linearly independent system of controls 
be complete in the subspace i.e. the relation holds, 

where V is a closure of the linear span (in the relevant norm). By (3.9), the 
system of waves 


u 


«;=n^.(.,T) = PF^/« 


turns out to be complete in Hi, i.e., one has = Hi- 

If T is such that \ 7 ^ 0, i.e., the waves moving from a do not cover 

the whole fl, then the control operator is injective [ 1 ] (in particular, this 
holds for T < To). In this case, W'^ preserves the linear independence, and 
turns out to be a linearly independent complete system in Hi- 


Convention 3. By this, we deal with this case and say to be a wave 

basis in the subspace Hi- Also, everywhere, system producing the 

wave basis, is chosen so that all /| G Domi?^. 
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As a consequence, the Gramm marices 

at - {(h.fGqS-i. n = i,2,... 

are nonsingular and invertible, whereas their entries can be represented via 
the controls: 

(at)a = . (3.10) 


• In the BCM, wave bases are used for hnding the projections of functions 
on the domains hlled with waves. 

Fix a positive ^ < T. Let Pj be the (orthogonal) projector in "Hj onto 
Such a projector cuts off functions: 

y in Gf 

0 in \ ni 

As an element of the subspace "Hi, this projection can be represented via the 
wave basis: 

N 

P^y = lim Pf^ = lim (3.11) 

k=l 

where P^^ projects in Pj onto the span and the column of coef- 

hcients {c^ ='■ is determined via the column {(?/, =: 

through the Gramm matrix by 

G = 

The limit is understood in the sense of the norm convergence in P^. 




3.3 Dual system 

• Denote Pj := {{x,t) \ 0 < T{x,a) <t < T}. 

A dynamical system associated with the problem 


vtt — c^Av = 0 

in Pj 

(3.12) 

n|t=r = 0, Vt\t=T = y 

in DJ 

(3.13) 

V = 0 

on X [0, T] 

(3.14) 
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Figure 2: Domain iFj 


is called dual to system aj; by u = v'^{x,t) we denote its solution. Owing 
to the FDI, such a problem turns out to be well possed for any y G "Hj. Its 
peculiarity is that the Cauchy data are assigned to the final moment t = T, 
so that the problem is solved in reversed time. 

The solutions to the original and dual problems obey the duality relation: 
for any / G and y G "Hj, the equality 

{uf(-,T),y)^r = UX)ri (3-18) 


is valid [HIEIEIE]. 

• With the dual system one associates an observation operator 

J-J, 

o'^y := ■ 


Writing (3.15) in the form {W'^f,y)^j. = {f,0^y)j^, we get an operator 
equality Hence, the dehnition of C'^ implies 


. 


(3.16) 


3.4 Projections of harmonic functions 


• Assume that ?/ in (3.13) is harmonic: y = a E TfJ obeys Aa = 0 in and 
is continuously differentiable up to cr^ C dfl^. A simple integration by parts 
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in (3.15) leads to 


{a,u^{-,T))^T = {0^a,f)^T = 

= [ {T-t)[a{j){R^- au{j)fi-f,t)] dTdt (3.17) 


’x[0,T] 


(see [SI El El HD])- 

Assume that C is chosen in accordance with Convention 

1^ and produces the wave basis C "HI. Then, representation (3.11) 

takes the form 


N 


P^a = lim 

k=l 

where satisfies the linear system 




(3.18) 


(3.19) 


with the Green matrix 




N 

*J=1 


= <2 


1-1 


{Rl^JS^f!){j,t) {R^ n){^,t)drdt 


T fp 


N 


(3.20) 


Jcrx[0,2T] j ij=\ 

and the right-hand side 

Bn = {(«. ul)n}Li, (a, ^ 

JaTx[0,T] L 


dTdt. 


With regard to (3.3),(3.4), and Convention]^ the latter can be written i 
the form 


m 




(T - () a(7)(fif S^/«)(7,() - o„(7)/|(7,() 


'<tJ’x[0,T] 


dVdt 


N 


k=l 

(3.21) 


determined by R?J" and, thus, relevant for the further use. 
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• Fix a positive ^ < T. The operator 


pS _ pT _ p5 

^(jl. ■ ^cr 


is the projector in 1-L onto the subspace it cuts off functions on the 

subdomain \ f2|. 

Choose systems {f^}k=i and {f^}k=i, which are linearly independent and 
complete in a nd respectively. Applying the (bounded) observation 


operator to (3.18), with regard to f we obtain 

N 


N 


O^Pja = jim cl^C^fl , O^Pja = lim^ 4 ,nC^ fk • 


k=l 


N^oo ■ 


k=l 


Subtracting, we arrive at the representation 


N 


k=l 


0^pL(^ = I^nC^ fk - cInC^ fi 




(3.22) 


For the future application to the inverse problem, a crucial fact is that its 
right-hand side is determined by the response operator. Indeed, if R'^ is 
given, one can 

1. choose the complete linearly independent systems C Pj and 

{fk}T=i P then, co mpose the Gramm matrices Qjf, by (3.20) 
and columns Bjf, by (3.21) 




2. solving system (3.19) with respect to hnd the coefficients 


3. determine C'^ by (3.6), compose the sum in (3.22) and, extending N, 
pass to the limit. 


3.5 Amplitude formula 

In what follows, we deal with the regular case T <T„. 

• Fix a positive ^ < T; let ?/ be a smooth function in fl. Return to the 


dual system (3.12)-(3.14) and put 


Vt\.T. = =: y\_ 
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in Cauchy data (3.13). Such a is of two specihc features: 

(i) it vanishes in so that suppj/^ is separated from a by the c-distance 
Therefore, by the hniteness of the wave propagation speed, vanishes in 
the space-time domain {(x, t) G iCj | t > {T — ^) + t{x, a)} and, in particular, 
one has 

= 0 forT — ^<t<T. (3.23) 

{a) is discontinuous: generically, it has jumps at the equidistant surfaces 
and r|. In particular, at the points G fl Tf, the value {ampli¬ 

tude) of the jump is 


2/i(a;(7,^ + 0)) = 2/(x(7,0)- 


(3.24) 


• In hyperbolic equations theory, the well-known fact is that discontin¬ 
uous Cauchy data initiate discontinuous solutions, the discontinuities prop¬ 


agating along characteristics. In our case of the wave equation (3.12), the 


jumps of Vt\^_rp = y± induce the jumps of in . In particular, there is a 
jump on the characteristic surface {(x, t) G Kji \ t = T — ^ t{x, a)} includ¬ 
ing its smooth part := {{x,t) G iCj | x G (on Fig 2, contoured by 

{15,16,17,18,15}). The jumps of on and of on the cross-section 
n = {(7, T — ^) I 7 G cr} (the line (15,16}) can be found by standard 
geometrical optics devices. For the latter jump, a simple analysis provides 


A I 


(7,^)|*=rJ!o = -/^(7,0l/(^(7,0), 7e 


a 


(see, e.g., [2] , [13] , [9] ). By (3.23), we have ui(-^(7,t)|* ^ = 0 that leads to 


u^^(7,T-^-0) = /3{-f,0y{x{-f,0), 'J ea. 


(3.25) 


Comparing (3.24) with (3.25), one can recall the well-known physical prin¬ 
ciple: jumps propagate along rays (here, a ray is = {x{'y, s) | 0 < s < ^}) 
with the speed c, the ratio of the jump amplitudes at the input and output 
of the ray (here, at x{'y,^) and x(7,0)) depending on the ray spreading. 

• Recalling the dehnitions of images and observation operator, one can 


write (3.25) in the form 


{O^Pl^y)h,0 = y{%0, (7.?)ey 


(3.26) 
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It is the so-called amplitude formula (AF), which plays a central role in solv¬ 
ing inverse problems by the BCM mw- It represents the image of function 
in the form of collection of jumps, which pass through the medium, absorb 
information on the medium structure, and are detected by the external ob¬ 
server at the boundary. 


• Now, let j/ = a be a harmonic function. Combining (3.22) with (3.26), 
we arrive at the key relation 


5(7. a = 

(7,0 6 sy 



|(7.() 


(3.27) 


As was noted at the end of section [3^ to hnd its right-hand side, it suffices 
to know the response operator. In particular, since the coordinate functions 


are harmonic, applying (3.27) to a = tt*, i = 0,... ,n one can recover their 
images tt* via R^. 


4 Determination of speed 

4.1 Procedure 

To solve the inverse problem, we just summarize our considerations in the 
form of the following procedure. Recall that the role of the procedure input 
data is played by operator . 

Step 1. Fix a. f < T. Applying the procedure 1. — 3. described at the end 


of section 3.4, hnd the right-hand side of (3.27) for a = vr^, tt^, ..., vr”' and, 

3 (7,0 ^ 


thus, get the images 7 r*( 7 ,^) for 7 G cr. 

Step 2. Varying hnd vr* on Sj. Then, recover the map Sj 
(r(7,0 ^ by the hrst representation in (2.3). The image of the map is 
RJ, so that the ray tube is recovered in hi. 

Step 3. Diherentiating with respect to hnd c by the second representation 
in (2.3). The pairs {(^( 7 , .^), c(a;( 7 , .^)) | (7,0 ^ constitute the graph of 


c in B^. 

Thus, the sound speed in the tube is determined. The following is some 
comments and remarks. 

• For applications in geophysics, by the obvious reasons, it is desirable 
to minimize the part of the boundary, on which the external observer has 
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to implement measurements. As is seen from (3.21), our procedure requires 
observations not only on a but on r\cr, whereas the knowledge of the bound 
c* just enables the observer to restrict measurements on aj. In principle, 
one can avoid the observations on r\cr by the use of the artihcial coordinates 
instead of the Cartesian tt*. Namely, one can choose the harmonic functions 
a^,..., a” obeying = 0, which separate points of the tube Bj at least 

locally. By this choice, in (3.21) one gets /^Txiqt] ~ /frx[or]- Therefore, 
possessing the values of f on (Tx[0,2T] (but not on the whole 0^^!), 
one can recover the images a® via the amplitude formula and use them for 
identifying the points of Sj in Thereafter, one recovers cl j.. However, it 
is not clear, whether this plan can provide workable numerical algorithms. 


• As was mentioned in |T3 the procedure j5] (sections 7.6, 7.7) enables 
one to determine c\^j. from observations on crx [0,2T] only, and, thus, pro¬ 
vides the strongest uniqueness result. However, its numerical implementation 
in the case of the time-domain inverse data seems to be rather problematic. 


4.2 Numerical testing 

Preparation of tests 

• We take 


H := {(x^, x'^) e I < 0} , T := {(x^, 0) € | — oo < < cx)} , 

a := {x G r I — L < x^ < L} , 


and consider a few concrete examples of the density p = in H. 

• We choose an appropriate hnite system of controls fk supported on a x 

[0,T], the system being the same for all examples except Test 1(6) where 
results with another spatial basis are presented for comparison. It is well 
adapted to constructing the systems of delayed controls: for intermediate 
f = fi, the shifts /f'(-,t) = ~ {T — ^i)) are in use. This enables one to 

reduce considerably the computational resources. 

• At each of the examples, we solve numerically the forward problems 
(1.1)-(1.3) with the hnal moments t = T and t = 2T for the controls 
fp and JS^ff‘ respectively. These problems are solved by the use of a 
semi-discrete central-upwind third order accurate numerical scheme with 
WENO reconstruction suggested in d- As a result, we get the functions 
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Rl^ entering in 


«S‘ = and uR'-' = 

( |3^ and ( |3^ . 


Controls 

The BCM uses a system of boundary controls fi, / 2 , ..., which belong to 
the Sobolev class: 


{/,Gi/'(rx[ 0 ,T])|/,( 7 ,t)|,=o = 0 } 

and constitute a basis in L 2 (r x [0, T]). We construct such a system from the 
products of elements of spatial and temporal bases, fki'Jyt) = 'iprnit), 
k = I + mN^, where Z = 0 : — 1, m = 0 : Nt — 1, and the basis dimension 

is N = N^Nt. 

In the case of the half-plane, we can keep under control only a part of 
the boundary and thus have to use localized basis functions. The simplest 
and good choice is a conventional trigonometric basis reduced to the interval 
[— 1 , 1 ] by an exponential cutoff multiplier 7 ( 7 ) = 1/(1 -|-exp ( 7 /s)) with a 
cutoff scale s, so that 


Ml) = vil - 1 ) vMl - 1 ) cos 


TT 


1 . Z -|-1 ., , 

2 + L—J(T-l) 


(4.1) 


where [-J is the integer part. The spatial basis functions are shown in the 
left panel of Fig. 

The temporal basis is constructed from the shifts of a tent-like function. 


9{t) 


d 

A 


{ 

■ A' 

\ 

1 1 — exp 


j In 


CQsh [^It^] cosh [^] 
cosh' [^] 


(4.2) 


so that 'ipmi't) = 0{t — mA — 6), where A = T/W, d is a smoothing parameter 
(when d —)■ 0 the function 9(t) gets a triangular shape), and S is an offset to 
ensure a negligible value of d(0). Such a shift-invariant basis (shown in the 
right panel of Fig. considerably reduces computational resources needed 
for the BCM-reconstruction. 


Regularization 


• In the course of determination of c by the procedure Step 1-3, we us e the 
above-prepared data for computing the entries of and in (3.21) and 
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(3.20). Then the system (3.19) is solved for a = tt^, tt 
are calculated by standard LAPACK routines 


° and the solutions 


, vr‘ 


Solving system (3.19), we have to apply a regularization procedure since 
the condition number of the Green matrix rapidly grows as its size N in¬ 
creases, see Figure]^ Because of unavoidable errors in matrix elements and 
right hand sides, the expansion coefficients also contain errors amplihed 
by ill-conditioned matrix. We use Tikhonov’s regularization to reduce fake 
oscillations caused by errors in expansion coefficients. The value of regular¬ 
ization parameter is selected to satisfy a desired tolerance for residual of the 
linear system. 

• One more operation, which produces unavoidable errors, is computation 


of the double limit in (3.27). The origin of the errors is the following. 


In (3.18), the projection Pja is a piece-wise smooth function in G, which 


has a jump at the surface T^. Therefore, the convergence of the sums in the 
right hand side not uniform near r|, and the Gibbs oscillations do occur in 
the summation process. These oscillations are transferred to the amplitude 


formula (3.27) and considerably complicate the determination of jump at 
t = T — whereas this determination is a crucial point of the algorithm. 

To damp this negative effect we apply the following procedure. The 
basis functions have hnite resolution of the order of spatial-temporal scales 
of the highest harmonic. All scales below the minimum ones are unreachable. 


therefore we average the result of expansions (3.27) over that minimum scales 
by convolution with some kernel 


+ O0 


+00 


{s)(7.*)= / d(' / 


(4.3) 


In our implementation the kernel — — is a product of conventional 

Gaussian kernels both for spatial and temporal variables. Such a procedure 
efficiently removes the Gibbs oscillations and, in fact, accelerates convergence 
of the expansions. The values of standard deviations in the Gaussian kernels 
should match the minimum spatial and temporal scales of the boundary 
controls to smooth out the oscillations. 


At the hnal step, the speed c is found by (2.3) with the help of numerical 


differentiation by the central hnite difference formula. 
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Numerical results 


Test 1. Let 

p(x\x^) = 1 + agi{x^)g 2 {x‘^), gk{x'") = exp 


{x^ — x^y 
2A^ 


(4.4) 


where a = 1, = 0, = —0.5, Ai = 0.5, A 2 = 0.5. The sound speed 

c = p~^ is shown on Figure [^together with exact semigeodesic coordinates 
and wave front at t = T = 1. We test the recovering procedure for two rather 


different spatial bases (the temporal basis (4.2) consisting of 16 functions is 
the same in both cases). 

a) In this subcase, we use spatial basis composed from localized trigonometric 


functions (4.1). A typical image of harmonic function x^ is shown in Figure 


1^ where we observe the Gibbs oscillations on the left plot and the smoothing 
effect of convolution (4.3) on the right one. The condition number of matrix 


(3.20) for .^ = T is 1.5T0^ and parameter of Tikhonov regularization for all 
linear systems is 1-10“^. The standard deviations of Gaussian kernels in (4.3) 


for ( 7 ,t) are = 0.1875 and at = 0. 

The mapping x( 7 , ^) is shown in Figure]^ The reconstruction error grows 
towards the ends of the localization interval 7 G (— 1 , 1 ) and for large values 
of ^ ^ T. 

The end result of the BGM is the sound speed recovered in the Gartesian 
coordinates. It is shown in central panel of Figure Relative errors of recon¬ 
struction in percents are shown in left panel of Figure]^ As is seen, although 
the reconstruction error quickly grows towards the ends of the localization 
interval and for large values of ~ T, in the most part of the domain covered 
by the direct rays from the boundary, the relative error does not exceed a 
few percents. 

b) Here we use a spatial basis composed from smooth tent-like functions 
as in (4.2). The condition number of matrix (3.20) for ^ = T is 5.9T0^ and 
parameter of Tikhonov regularization for all linear systems is 1-10“®. The 


standard deviations of Gaussian kernels in (4.3) for ( 7 ,t) are a^ = 0.125 and 
at = 0.0625. The recovered speed of sound is shown in right panel of Figure 
while its relative errors are shown in right panel of Figure 
We may conclude that both of these bases provide similar quality of 
reconstruction of the order of several percents in most part of the domain. 
The advantages of the tent-like basis are smaller condition number of the 
system matrix and the same spatial scale of all basis functions. The effect 
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of lower accuracy of reconstruction along the lateral boundaries in the case 
(b) is due to narrower support (smaller value of L) of the boundary controls 
compared to the case (a). 

Test 2. For the second test, we take 


p{x^, = 1 — O.Sx^ + 0.0625 ~ ® 


dg2^ 

dx"^ 


where a = 0.25, x^ = 0, = —0.5, Ai = 0.5, A 2 = 0.25. The corresponding 

sound speed has a background value 1 and two variations of the order 30% 
of its boundary value. 

We use T = 1.5 and the basis with 16 spatial (trigonometric) and 32 


temporal functions. The condition number of matrix (3.20) is shown in Fig- 
it grows as This is a consequence of the strong ill-posedness of 


ure 


the inverse problem under consideration, and such a growth constrains the 
maximal depth of reconstruction (determined by errors in right hand sides 


(3.21)), which is possible for the given part of the boundary. In computa¬ 
tions, the parameter of Tikhonov regularization for all linear systems is hxed 
and equal to 1-10“^. The standard deviations of Gaussian kernels in (4.3) 
for (7,f) are = 0.1875 and at = 4.6875-10“^. For ^ ^T, the error in the 
expansion coefficients a| grows up and we had to increase a^ to the value 0.5 
for smoothing out the large scale fake oscillations from low spatial harmonics. 
Such an over-smoothing reduces the accuracy of the recovering for ^ ^ T. 

The recovered speed of sound in the Cartesian coordinates is shown in 
and relative errors of reconstruction in percents are shown in 
Thus, in the most part of the domain Qcr covered by the direct 


Figure 10 


Figure 11 


rays coming from a, the relative error does not exceed a few percents. 


Test 3. Here we take 


p{x^,x‘^) = 1 — 0.5x^ -f 0.0625 {x^Y + agi{x^) (l — x^) 


dg2{x^) 
dx^ ' 


where a = 0.25, x^ = 0, x^ = —0.5, Ai = 0.5, A 2 = 0.25. In contrast to the 
case 2, the corresponding speed of sound has rather strong variations (the 
ratio of maximum to minimum value is about 2.5). 

Again, we take T = 1.5 and use the basis with 16 spatial and 32 temporal 
functions. The condition number of matrix (3.20) for .^ = T is 2-10^ and in 


calculations the parameter of Tikhonov regularization for all linear systems 
is 1-10“^. The standard deviations of Gaussian kernels in (4.3) for ( 7 ,t) are 
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= 0.1875 and at = 4.6875-10“^. Again, for ^ T we had to increase a^ 
to valne 0.625 to smooth ont large scale fake oscillations from the low spatial 
harmonics. 

The recovered speed of sonnd in the Cartesian coordinates is shown in 


Fignre 12, and relative errors of reconstrnction in percents are shown in 
Fignre[T^ 

Test 4. To test the recovering algorithm for the case of sonnd speed qnickly 
varying along the bonndary we set 


p{x^,x^) = 1 - ag2{z^) 


dx^ 


= cos(0)x^ + sin(0)(a:^ + 0.25), 

= — sm{(j))x^ + cos{(j)){x‘^ + 0.25), 

where a = 0.25, x^ = 0, = 0, Ai = 0.375, A 2 = 0.25, (j) = 7r/12. 

We take T = 1 and nse the basis with 16 spatial and 32 temporal fnnc- 
tions. The condition nnmber of matrix (3.20) for .^ = T is 1.28T0^ and in 


calcnlations the parameter of Tikhonov regnlarization for all linear systems 
is 1-10“^. The standard deviations of Ganssian kernels in (4.3) for ( 7 ,t) are 


a-y = 0.125 and at = 3.125-10“^. Again, for ^ T we had to increase a^ to 
the valne 0.25 in order to decrease large scale fake oscillations from the low 
spatial harmonics. 


The recovered speed in the Cartesian coordinates is shown in Figure 14 


and relative errors of reconstruction in percents are shown in Figure 15 


We observe rather large oscillations in the speed values emerging at = 
—0.6, and to demonstrate the origin of these oscillations we also show the 
results of a pseudo-reconstruction. The latter means the use of a conventional 
recovering procedure, in which all the matrix elements (products (M-^ynb)) 
are computed via the solutions found by solving the forward problem with 
the given (known) speed prohle. Such products are much more accurate than 
the ones found via the inverse data, and therefore the errors in the expansion 
coefficients in (3.27) are greatly reduced in the pseudo-reconstruction. This 


leads to much better quality of recovering far from the boundary and clearly 
shows the effect of ill-posedness of the reconstruction procedure. 

Test 5. The purpose of the last test is to check the ability of BCM to work 
with strong gradients in the recovered quantities. We prepare the density 
of medium as a slightly smoothed wedge with density p = 5 included in the 


homogeneous background with the constant density p = 1, see Figure 16 


24 












We take T = 1 and use the basis with 16 spatial and 32 temporal func¬ 
tions. The condition number of matrix (3.20) for .^ = T is 2.07T0^ and in 
calculations the parameter of Tikhonov regularization for all linear systems 
is 1.0T0“^. The standard deviations of the Gaussian kernels in (4.3) for ( 7 , t) 
are = 0.125 and at = 3.125-10“^. Again, for^^T we had to increase a^ 
to the value 0.325 in order to decrease large scale fake oscillations from the 
low spatial harmonics. 

The recovered density in the Cartesian coordinates is shown in Figure 
and relative errors of reconstruction for c{x) in percents are shown in 
The location of the wedge is recovered with good accuracy and 
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Figure 17 


without systematic shifts, the maximum of recovered density is 4.68 that 
corresponds to the relative error about 7%. The smearing of discontinuities 
is quite reasonable taking into account the spatial scales of boundary controls 
and the smoothing ranges, whereas the errors at the discontinuities are big, 
as it has to be. We may suggest that if the resolution of boundary controls 
is not enough to resolve spatial scales of the medium inhomogeneities, then 
the method will recover a smoothed averaged prohle, which can be further 
improved by another high resolution methods. 


Comments 

• Our results on the numerical speed determination from the time domain 
data given at a part of the boundary demonstrate that the BCM-algorithm 
is workable and provides good reconstruction in the domain covered by the 
normal acoustic rays. 

• The key step of the algorithm is inversion of a big-size Gram matrix 
(AT 10 ^ — 10 ^), which consists of the inner products of waves initiated 
by rich enough system of boundary controls. As is typical in multidimen¬ 
sional (strongly ill-posed) inverse problems, the condition number of this 
matrix rapidly grows as one extends the number of controls and/or the ob¬ 
servation time. Therefore, to increase the depth of reconstruction one has to 
use controls acting from a larger part of the boundary or decrease errors in 
the input inverse data. 

• The number and shape of boundary controls determine the spatial reso¬ 
lution of the reconstruction procedure. The BGM is able to work with low 
number of spatial controls: in such a case it provides an ‘averaged’ prohle. 
As we hope, such a prohle can be used as a starting approximation for high 
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resolution iterative reconstruction methods B 

In the future work, we plan to evaluate the influence of external noises 
in the inverse data on the quality of the BCM reconstruction, and apply the 
method to another domains and more realistic sound speed prohles. 
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Figure 3: Basis of boundary controls: spatial functions 0 ^( 7 ) with Z = 0 : 16 
and s = 1/32 (left) and temporal functions with m = 0 : 15 and T = 1, 
Nt = IQ, d = A/QA (right). 
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Figure 4: Test 1. Speed of sound c{x) (color), exact semigeodesic coordinates 
(mesh), and wave front at t = T (line). 
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Figure 5: Test 1. Image vr^( 7 ,.^): expression (3.27) (left) and its smoothed 
version (4.3) (right). 



Figure 6: Test 1. Reconstructed mapping x = xi^j,^), the exact values are 
shown by black mesh. 



Figure 7: Test 1. Speed of sound c{x) in the domain hlled by waves initiated 
from a: left - exact values, central - recovered values with trigonometric 
spatial basis, right - recovered values with tent-like spatial basis. 


30 



























Figure 8: Test 1. Map of relative errors (in percents) of the recovered sound 
speed c{x): left - with trigonometric spatial basis, right - with tent-like spatial 
basis. 


Figure 9: Test 2. The condition number of matrix (3.20) of scalar products 
as function of the probing time 
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Figure 10: Test 2. Speed of sound c(x) in the domain filled by waves initiated 
from a: left - exact values, right - recovered values. 



Figure 11: Test 2. Map of relative errors of the recovered sound speed in 
percents. 
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Figure 12: Test 3. Speed of sound c(x) in the domain filled by waves initiated 
from a: left - exact values, right - recovered values. 



Figure 13: Test 3. Map of relative errors of the recovered sound speed in 
percents. 
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Figure 14: Test 4. Speed of sound c(x) in the domain filled by waves initiated 
from a: left - exact values, central - recovered values, right - pseudo-recovered 
values. 



Figure 15: Test 4. Map of relative errors (in percents) of the recovered sound 
speed c(x): left - usual reconstruction, right - pseudo-reconstruction. 
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Figure 16: Test 
initiated from a: 


5. Density of medium p{x) in the domain filled by waves 
left - exact values, right - recovered values. 



Figure 17: Test 5. Map of relative errors of the recovered sound speed in 
percents. 
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